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Perturbative expansions of several small Wilson loops are computed through next-to-next-to- 
leading order in unquenched lattice QCD, from Monte Carlo simulations at weak couplings. This 
approach provides a much simpler alternative to conventional diagrammatic perturbation theory, 
and is applied here for the first time to full QCD. Two different sets of lattice actions are considered: 
one set uses the unimproved plaquette gluon action together with the unimproved staggered-quark 
action; the other set uses the one-loop-improved Symanzik gauge-field action together with the 
so-called “asqtad” improved-staggered quark action. Simulations are also done with different num¬ 
bers of dynamical fermions. An extensive study of the systematic uncertainties is presented, which 
demonstrates that the small third-order perturbative component of the observables can be reliably 
extracted from simulation data. We also investigate the use of the rational hybrid Monte Carlo 
algorithm for unquenched simulations with unimproved-staggered fermions. Our results are in ex¬ 
cellent agreement with diagrammatic perturbation theory, and provide an important cross-check of 
the perturbation theory input to a recent determination of the strong coupling «Mg {Mz) by the 
HPQCD collaboration. 


I. INTRODUCTION 

A key ingredient in many high-precision applications 
of lattice QCD is the use of perturbation theory, in order 
to match lattice discretizations of actions and observ¬ 
ables to their counterparts in continuum QCD. An im¬ 
portant example is the determination of the strong cou¬ 
pling ctyig(Mz), where perturbative expansions of short- 
distance quantities, such as small Wilson loops, can be 
used to extract the value of the coupling from simula¬ 
tion measurements [I, 0. Lattice perturbation theory by 
Feynman diagram analysis is extremely difficult however, 
because the lattice regulator results in Feynman rules 
that are exceedingly complicated; moreover, there is a 
proliferation of diagrams that are not present in the con¬ 
tinuum. A further challenge is that these perturbative- 
matching calculations must generally be carried out at 
next-to-next-to-leading order (“NNLO,” which is gener¬ 
ally equivalent to two-loop Feynman diagrams), if one is 
to obtain results of a few-percent precision |j. 

Although many parts of diagrammatic lattice pertur¬ 
bation theory have been automated with the help of com¬ 
puter codes 000 , higher-order calculations remain 
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very challenging, and very few NNLO lattice calculations 
have been done (see, e.g., Refs. [2, 000) The algebraic 
burden is particularly heavy for the highly-improved ac¬ 
tions that are now commonly used in numerical simu¬ 
lations. A particularly important case is the tree-level 
0(a 2 )-improved staggered-quark action (where a is the 
lattice spacing) 0, together with the 0(o 2 V-accurate and 
one-loop Symanzik-improved gluon action [5[,|To|], both of 
which are tadpole improved (and which are hereafter col¬ 
lectively referred to as the “asqtad” actions). The “asq¬ 
tad” actions are currently being used by the MILC col¬ 
laboration to generate unquenched ensembles with three- 
flavors of sea quarks EH , which have been used by several 
groups for a wide variety of physics applications, includ- 
ingq uantities requiring higher-order perturbation theory 
000 

Given the central role of perturbative matching in 
lattice QCD, alternatives to diagrammatic perturbation 
theory are desirable. One approach [3 jl3j] is to use 
Monte Carlo simulations at weak couplings, where the 
theory enters the perturbative phase (at finite volume). 
Simulation measurements of a particular observable are 
done at several values of the coupling, and the result¬ 
ing data are fit to an expansion in the coupling; the fit 
yields numerical values for the perturbative coefficients, 
without Feynman diagrams, and with little or no analytic 
input. 

This Monte Carlo approach has been successful in com- 
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puting the perturbation series of a number of quantities 
in pure gauge theories [H HI HI HI In particular, 
perturbative coefficients for a number of small Wilson 
loops for the plaquette gluon actionjwere computed to 
NNLO using this method in Ref. 0 , and the results 
were subsequently reproduced by diagrammatic pertur¬ 
bation theory (see Ref. Jj). 

In this paper we make the first applications of this 
Monte Carlo method to lattice actions with dynamical 
fermions 0 . The perturbation series of a number of 
small Wilson loops are computed through NNLO, for 
two different lattice QCD actions with dynamical stag¬ 
gered quarks: the Wilson plaquette gluon action with the 
unimproved-staggered quark action, and the “asqtad” ac¬ 
tions. The latter calculation is of particular relevance be¬ 
cause it provides a consistency check of the NNLO Wilson 
loop expansions used in the determination of (Mz) 
by the HPQCD collaboration, in Ref. 

We note that another simulation approach to lattice 
perturbation theory has also been developed, based on an 
explicit perturbative expansion of the simulation equa¬ 
tions themselves El, and has recently been applied to 
(unimproved) lattice actions with dynamical fermions 

The Monte Carlo method considered here has the ad¬ 
vantage that it can be done using conventional simulation 
codes, although it is very advantageous to a dop t twisted 
boundary conditions (TBC) [[|. 0, f20, 0, [22| in order 
to eliminate zero modes, and to suppress nonperturba- 
tive finite-volume effects due to transitions between Z( 3) 
phases 0 ; fortunately, TBC can be implemented in ex¬ 
isting simulation codes with relatively little effort. 

This simulation approach to unquenched perturbation 
theory is also very efficient, since the simulations can be 
done on very small lattices (8 4 volumes are used here), 
thanks to the use of TBC. The simulations are also done 
here for different numbers rif of dynamical fermions, 
which explicitly demonstrates that the data are sensitive 
to the n/-dependence of the perturbative coefficients. 

Although small Wilson loops are the simplest observ¬ 
ables for perturbative analysis by Monte Carlo simula¬ 
tions, because they are such ultraviolet quantities (as 
well as being gauge invariant), the studies presented here 
nonetheless provide important benchmarks for future ap¬ 
plications to other quantities; for instance, similar meth¬ 
ods were used to analyze the static-quark propagator in 
pure-gauge backgrounds in Ref. EH, and the Fermilab 
heavy-quark propagator was analyzed in Ref. El- 

While the expansion coefficients can be extracted with 
far less effort from Monte Carlo simulations, care must 
be taken to control all of the systematic uncertainties, 
in order to reliably extract the small part of the signal 
corresponding to the higher-order terms in the perturba¬ 
tive expansion. It is also crucial to use an expansion in 
a renormalized coupling, rather than in the bare lattice 
coupling ai a t = < 7 2 /(47t), for which perturbation theory 
is very poorly convergent 0|. We use an expansion in 
the coupling ay(q*) defined by the static potential, along 


with an estimate of the optimal scale q* for a given quan¬ 
tity 0 , [24j. Although one can in principle define the 
renormalized coupling ay at a given bare lattice coupling 
cq a t from simulation measurements of the static poten¬ 
tial, we rely here instead on existing NNLO determina¬ 
tions of the relationship from diagrammatic perturbation 
theory [?). 

We did simulations at couplings ay < 0.1, hence sta¬ 
tistical and systematic errors must be much less than 
max(ay) ~ 10~ 3 , for determination of the third-order 
coefficients. The ensemble sizes were chosen in order to 
reduce the statistical errors to the desired level. There 
are four major sources of systematic error in the simula¬ 
tion algorithms and in the data analysis: i) transitions 
between the Z( 3) center phases of the gauge field action; 
ii) finite step-size errors in the simulation equations; iii) 
other algorithmic systematics, such as the precision of the 
matrix inversion; and iv) fitting and truncation errors in 
the perturbative expansion. 

In this study simulations were done with TBC to sup¬ 
press transition between the center phases, as indicated 
above. To eliminate the finite step-size errors we pursued 
one of two strategies: in the case of the unimproved ac¬ 
tions, we used the rational hybrid Monte Carlo algorithm 
(RHMC) pH 12(1 : for the “asqtad” actions we used the 
R algorithm )27l| . and did simulations at several differ¬ 
ent molecular-dynamics step sizes, the results of which 
were extrapolated to zero step size. To extract the per¬ 
turbative coefficients the data were analyzed using con¬ 
strained curve fitting techniques l2Sj . which provide an 
elegant procedure for incorporating our a priori expecta¬ 
tion that the perturbative expansion is well behaved, and 
which readily allow the maximum amount of information 
to be extracted from the data. 

The rest of the paper is organized as follows. The ac¬ 
tions and simulation parameters are detailed in Sect. nn 
along with the perturbation theory input which we use 
to extract the renormalized coupling ay at a given bare 
lattice coupling, from the simulation data. A detailed 
analysis of the systematic uncertainties is presented in 
Sect. If 111 including the use of TBC to control finite- 
volume effects. Results for the two actions are reported 
in Sect. inn and are compared with NNLO diagrammatic 
perturbation theory. Conclusions and prospects for fur¬ 
ther work are presented in Sect. El 


II. ACTIONS AND SIMULATION 
PARAMETERS 

A. Perturbative expansions 

We have done unquenched simulations for two sets of 
actions: the first set used the unimproved Wilson pla¬ 
quette action with unimproved-staggered fermions, while 
the second set used 0{a 2 )-improved gluon and staggered- 
quark actions. We provide the simulation parameters and 
perturbation theory input for each set of actions in the 
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two following subsections. 

A basic input parameter for the simulations is the bare 
lattice coupling oq at = g 2 /{4n), related to the usual sim¬ 
ulation parameter (3 in the case of the Wilson plaquette 
action, for example, according to (3 = 6 /g 2 = 6/(47rcq at ). 
However the bare coupling is not suitable for perturba¬ 
tive expansions 12.'ll and we use instead a renormalized 
coupling avia), defined by the static potential according 

to |E El 

s 44.MS) (1) 

3 q z 

A (truncated) perturbative expansion for the logarithm 
of a Wilson loop of size R x T is used, owing to the 
perturbative perimeter law, 

1 N 

~ r '(R I T) ~ y^ c rt;i?TQ^(gfl T ), (2) 

' n= 1 

where N is the order at which the series is truncated. 
We aim to measure the coefficients through third-order, 
hence fits must be done with N > 3. The characteris¬ 
tic scale q* RT for each observable is determined according 
to the BLM method (which estimates the typical mo¬ 
mentum carried by a gluon in leading-order diagrams) 

B0J- 

The connection between the simulation input «i a t and 
the ay coupling is therefore required. Although it is 
possible, in principle, to extract this relation by directly 
measuring the static-quark potential in Monte Carlo sim¬ 
ulations, the connection has already been computed in 
diagrammatic perturbation theory through NNLO for a 


variety of actions, including those we consider here, in 
Ref. [Z]. We use the results of Ref. 0] to provide the 
three leading orders of the expansion for the lxl pla¬ 
quette, as well as the scales q* RT for all the Wilson loops 
(the later only requires a relatively straightforward one- 
loop calculation). 

The diagrammatic input provided by the expansion 
of the lxl loop simplifies our analysis, and still al¬ 
lows for highly-nontrivial applications of the Monte Carlo 
method; here we demonstrate the utility of the method 
by computing the NNLO perturbative expansions of sev¬ 
eral larger Wilson loops, and thereby also provide a 
valuable consistency check of the NNLO diagrammatic 
calculations which were used to obtain a^jg (Mz) in 
Refs. 0,0. 

Given this input from diagrammatic perturbation the¬ 
ory, the Monte Carlo method proceeds as follows. Simu¬ 
lations are done at several small values of the bare cou¬ 
pling aiat (at which the lattice theory is in the pertur¬ 
bative phase at a given finite volume). For each bare 
coupling one measures the average plaquette (W'h)mc, 
whose perturbative expansion is assumed, and the quan¬ 
tities of interest, whose perturbative expansions are to be 
determined; in our case these latter quantities are larger 
Wilson loops, (Wrt) mc- The numerical value of the 
renormalized coupling ay at the scale q*± is obtained 
by substituting the measured value (Wh)mc into its 
known third-order expansion. Once the numerical value 
of av{qh) has been determined, the couplings at the 
other relevant scales q RT can be computed using the uni¬ 
versal second-order beta function, plus the known third- 
order correction in the V scheme according to 


a v {q) 


47r 

Po ln(<7 2 /Ay) 


f3i In [ln(g 2 /Ay)] ft 

Po Hq 2 / a v) Po ln2 (<l 2 / A v) 



P 2 VP 0 

Pt 



( 3 ) 


where av{<lii) is traded for the intrinsic scale Ay at 
the given bare coupling. The beta-function coefficients 
are (3q = 11 — |n/, (3\ = 102 — ^n/, anc ^ P' 2V = 
4224.18 — 746.006n/ + 20.8719n.j. Fits to the expansion 
Ecp 0 then yield the perturbative coefficients for the 
larger Wilson loops; details on the fitting procedure are 
described in Sect. mm 

B. Unimproved actions 

The first set of actions we consider are the Wilson pla¬ 
quette gluon action, 

C’=»E( 1 -w w 


where P Rl/ is the lxl plaquette and (3 = 6 /g 2 , together 
with the unimproved staggered-quark action 

‘-’stagg — +m 0 | x(x), (5) 

x \ n / 

where is the standard lattice derivative, and = 

(—1) X1+ -- +X ''- 1 is the usual staggered-quark phase. 

Two sets of simulations were done for the unimproved 
actions; the simulation parameters are summarized in 
Table d One set was done for a single dynamical fla¬ 
vor, n/ = 1, and another was done for three degener¬ 
ate dynamical flavors, n/ = 3. In each case, simula¬ 
tions were done for seven values of the bare coupling, 
with bare-quark mass m 0 a = 0.2, on 8 4 lattices with 
twisted boundary conditions (the boundary conditions 
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Number of flavors: nf = 1 

P 

«lat = 6/47T/3 

(Wn) 

Measurements 

Acc. Rate 

OLviqh) 

aAy 

11.0 

0.04341 

0.805569(41) 

849 

92% 

0.05574 

2.1 x 10" 8 

13.5 

0.03537 

0.843964(34) 

902 

90% 

0.04292 

9.1 x 10 -6 

16.0 

0.02984 

0.869561(26) 

992 

87% 

0.03496 

4.0 x 10 -7 

19.0 

0.02513 

0.890999(21) 

981 

81% 

0.02860 

9.2 x 10" 9 

24.0 

0.01989 

0.914413(16) 

1066 

69% 

0.02198 

1.7 x 10 -11 

32.0 

0.01492 

0.936259(11) 

1180 

58% 

0.01605 

7.2 x 10“ 16 

47.0 

0.01016 

0.956886(7) 

1261 

53% 

0.01067 

4.2 x 10~ 24 


Number of flavors: rif = 3 

P 

«iat = 6/47T/3 

(Wn) 

Measurements 

Acc. Rate 

av(<7n) 

aAy 

11.0 

0.04341 

0.807064(36) 

821 

89% 

0.05551 

3.9 x 10 _b 

13.5 

0.03537 

0.844847(26) 

641 

86% 

0.04283 

1.1 x 10" 6 

16.0 

0.02984 

0.870205(21) 

740 

84% 

0.03490 

2.9 x 10 -8 

19.0 

0.02513 

0.891387(18) 

738 

77% 

0.02859 

3.8 x 10” 10 

24.0 

0.01989 

0.914647(13) 

791 

66% 

0.02197 

2.7 x 10 -13 

32.0 

0.01492 

0.936402(10) 

840 

56% 

0.01605 

2.5 x 10” 18 

47.0 

0.01016 

0.956950(6) 

902 

52% 

0.01066 

8.5 x 10 -28 


TABLE I: Simulation parameters for the unimproved plaquette gluon action with the unimproved staggered-quark action. Two 
different sets of simulations were done for the indicated number of flavors; in both cases 8 4 lattices with twisted boundary 
conditions were used, with bare-quark mass mod = 0.2. The RHMC algorithm was used with step size At = 0.01, with 50 
molecular-dynamics steps per trajectory; the table shows the acceptance rate for the accept/reject step at the end of each 
trajectory. 


Loop 

Qrt 

Loop 

* 

Qrt 

1 x 1 

3.40 

2x2 

2.65 

1x2 

3.07 

2x3 

2.56 

1x3 

3.01 

3x3 

2.46 


TABLE II: The optimal momentum scales q^ T for selected 
RxT Wilson loops for the unimproved actions Q- 


are discussed in Sect. IIII All . The configurations for these 
actions were generated using the RHMC algorithm, de¬ 
scribed in more detail in Sect. ITTT"H1 we used a time-step 
At = 0.01 with 50 molecular-dynamics steps between ac¬ 
cept/reject tests. Based on an autocorrelation analysis, 
described in Sect. rrrm 40 trajectories were skipped be¬ 
tween measurements. 

For the unimproved actions the expansion coefficients 
of the average plaquette to NNLO are given by UM 

c 1;11 = 1.04720(0), 

c 2; ii = -1.2467(2) - 0.06981(5)n/, 

c 3;11 = -1.778(7) + 0.464(27)71/ + 0.00485(0)7^, (6) 

where the uncertainties in the coefficients are statisti¬ 
cal errors which arise from numerical integration of the 
multi-loop Feynman diagrams using a Monte Carlo tech¬ 
nique (the error of “0” in ci ; n indicates that the integra¬ 
tion error is in the sixth digit in that case). The relevant 
scales for the Wilson loops are given in Table UH 


C. Improved actions 

The one-loop Symanzik-improved gluon action [M3 
we use follows Refs. ITU iTT] for tadpole improvement 

SgZ = /3 P i E d - p ^)+ Ed- R ^) 

+ /^pg El (1 — C/b±A±a-)> (7) 

x;n<i'<<j 

where R is the 1x2 rectangle, is the lxlxl 
“corner cube” (see Ref. [01), and where the couplings 
are given by 

/3 P i = =-dE(l + 0.4805a s ), 

g ZUtig 

/3 pg = -%0.03325a s , (8) 

u 0 

with a s here defined by the (first-order accurate) expres¬ 
sion 

a s = —41n(ito)/3.0684. (9) 

The one-loop couplings in S'Z correspond to the 
average-plaquette definition of the gluon mean field 

u 0 = (TFn) 1/4 . (10) 

The tree-level 0(a 2 )-accurate improved staggered- 

quark action we simulate was derived in Ref. 9;j, and 
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p 

Input Mo 

Measured mo 

C^lat 

<M<7iT) 

aAv 

Measurements 

9.5 

0.91690 

0.916922(98) 

0.08377 

0.12709 

5.8 x 10~ 2 

461,457,836,1335 

11.0 

0.93166 

0.931687(73) 

0.07234 

0.10113 

2.0 x 10” 2 

281,633,1090,1038 

13.5 

0.94704 

0.946986(55) 

0.05895 

0.07604 

3.2 x 10~ 3 

290,635,1122,1078 

16.0 

0.95661 

0.956614(53) 

0.04974 

0.06106 

5.0 x 10“ 4 

296,634,1145,1100 

19.0 

0.96433 

0.964311(38) 

0.04188 

0.04951 

5.5 x 10' 5 

298,643,1123,1172 

24.0 

0.97243 

0.972453(29) 

0.03316 

0.03765 

1.3 x 10 -6 

308,645,1191,1113 

32.0 

0.97978 

0.979785(23) 

0.02487 

0.02727 

3.3 x 10' 9 

315,652,1205,1204 

47.0 

0.98652 

0.986526(15) 

0.01693 

0.01797 

3.8 x 10~ 14 

215,709,1227,963 

80.0 

0.99220 

0.992208(13) 

0.00995 

0.01029 

5.3 x 10~ 25 

472,489,974,1527 


TABLE III: Simulation parameters for the “asqtad” actions on 8 4 lattices with twisted boundary conditions. The number of 
flavors is n/ = 1, with bare-quark mass mo a = 0.1. Simulations at each bare lattice coupling were done at four values of the 
-R-algorithm step-size, At = 0.005,0.01,0.02, and 0.03. The measured uo and avilii) are only shown for the ensembles with 
At = 0.005. The number of measurements at each of the four step sizes is given. 


Loop 

Qrt 

Loop 

* 

Qrt 

1 x 1 

3.33 

2x2 

2.58 

1x2 

3.00 

2x3 

2.48 

1x3 

2.93 

3x3 

2.37 


TABLE IV: The optimal momentum scales q%r for selected 
RxT Wilson loops for the “asqtad” actions Q. 


is also used in the three-flavor simulations by the MILC 
collaboration m 


surements. 

For the “asqtad” actions the expansion coefficients of 
the average plaquette to NNLO are given by 0, 3Cj] 

c 1;11 = 0.76710(0), 

c 2;11 = -0.5945(2) - 0.07391(2)71/, 

c 3;11 = -0.589(38) + 0.600(2)71/ + 0.00774(0)r^, (12) 

The relevant scales for the Wilson loops are given in Ta¬ 
ble nvi 




Z! ’&*(*) ( a m ~ 


■ mo 


x(z); 

( 11 ) 


see Ref. Q for the definition of the “smeared” derivative 
operator A^. Tadpole improvement of is defined 

by replacing each link U^ in the action by U^/uq, but 
only after adjacent pairs of identical links = I) 

are eliminated; the tadpole weights for the individual link 
paths in the quark action were also detailed by the MILC 
collaboration Ha- 

Simulations were done for the “asqtad” actions only 
for a single dynamical flavor, rif = 1. Ensembles were 
generated at nine values of the bare coupling, on 8 4 lat¬ 
tices and with bare-quark mass mod = 0.1; simulation 
parameters are given in Table urn 

The R algorithm was used in these simulations, 
and ensembles were generated at each bare coupling 
at four different molecular-dynamics step sizes, At = 
0.005,0.01,0.02, and 0.03, with the number of steps per 
trajectory in the four cases given by 100, 50, 25, and 15, 
respectively. The measured Wilson loop data are extrap¬ 
olated to At = 0 at each bare coupling using constrained 
curve fitting (cf. Sect. mrm The simulation value of the 
tadpole factor uq is determined self-consistently at each 
bare coupling by iteration during the thermalization pro¬ 
cess; the “final” input value was verified for consistency 
with the final measured value of Mo, as shown in Ta¬ 
ble EH Based on an autocorrelation analysis, described 
in Sect. imn we skipped 20 trajectories between mea- 


III. SYSTEMATIC EFFECTS 

The various systematic effects which must be con¬ 
trolled in order to extract higher-order expansion coef¬ 
ficients were enumerated in the Introduction, and in the 
following four subsections we consider these in turn. 


A. Z( 3) phases and twisted boundary conditions 

The SU(3) co ior gauge-field action is invariant under the 
transformation 

—> zU/j,(x), Vx 9 x ■ jl = constant, (13) 

where z € [1, e* 27r / 3 , e i47r / 3 ] is an element of the center- 
subgroup Z( 3). Although all closed Wilson loops are 
invariant under this transformation, the Polyakov line in 
the /i direction is sensitive to the phase; moreover, the 
perturbative expansion of Wilson loops can be spoiled 
by the formation of domains of different center phases 
Ha . While the fermion action breaks the center sym¬ 
metry, the action “cost” for reaching different phases, at 
the lattice volume considered here, is too small to pre¬ 
vent frequent transitions between phases, at least when 
periodic boundary conditions (PBC) are used. This can 
be seen in Fig.E where a scatter plot and run-time his¬ 
tory of the “temporal” Polyakov line are shown for the 
unimproved actions with PBC, on a 4 4 lattice at (3 = 16, 
well into the deconfined phase of the theory. 
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FIG. 1: Scatter plot and run-time history of the “temporal” Polyakov loop Pt on a 4 4 lattice at /3 = 16, using periodic 
boundary conditions. Results are shown for the unimproved actions using the RHMC algorithm, with n/ = 1, and bare-quark 
mass moo = 0.2; 100 trajectories were skipped between measurements. The simulations were started with all links set to the 
identity, U ll (x) = I, and the history is shown starting from this initial configuration. 


Although the scatter plot shows that the action with 
dynamical fermions has a preference for configurations 
near the non-trivial phases z = e* 2 ’ 1 ’/ 3 and 0 = e l4,r//3 
(despite the fact that the simulation was started with all 
links set to the identity), owing to the breaking of the 
Z( 3) symmetry, transitions between the various phases 
occur frequently, which would prevent the extraction of 
perturbative expansions, except on very large lattices 

Fortunately, the Z{ 3) transitions can be largely elimi¬ 
nated by using twisted boundary conditions (TBC) |l3j, 
which can be easily implemented using existing simula¬ 
tion codes. In the case of the gauge fields, there is no 
change to the link variables inside the lattice, and link 
variables across the twisted lattice boundaries are com¬ 
puted according to 

U^x + Lu) = (14) 

where L is the lattice length, and the are a set of 
constant “twist” matrices which obey the algebra [fj 

= 2 n 3 oc/. (15) 


Twists must be applied across at least two boundaries, 
else the effect of the twist matrix can be removed by 
a field redefinition. We choose to impose twists across 
all three “spatial” lattice boundaries, with ordinary peri¬ 
odic boundary conditions across the “temporal” bound¬ 
ary m ■ The following explicit representation of the twist 
matrices was used in our simulations 


n.T — 


' 0 

1 

o' 


■ g-i27r/3 

0 

0 

0 

0 

1 

II 

a 

0 

1 

0 

1 

0 

0 

0 

0 

e +i27r/3 


= ntn v = 


o 
i 

o e 


o 

0 

+i27r/3 


e ~i2Tv/3 

o 

0 


(16) 


To apply twists to the unquenched theory the quark 
field fields x( x ) must become 3x3 color matrices l20li2ll : 
these additional fermionic degrees-of-freedom amount to 
a new set of three degenerate “flavors” (Parisi introduced 
the term “smells” |2(J, to distinguish these copies from 
physical quark flavors). To compute the fermion fields 
across the twisted lattice boundaries one then imposes 
the boundary conditions 


The gauge action and observables are therefore periodic 
with TBC, but with period 3L, and with the important 
additional benefit that zero modes are eliminated |{|. 


X{x + Lv) = e™ ,3 VL v x{x)tii, (17) 


where the phase e”’/ 3 makes the quark field anti-periodic 
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FIG. 2: Unimproved action scatter plot and run-time history, with the same simulation parameters as in Fig.0 but here using 
twisted boundary conditions (TBC) in the three “spatial” directions. 


(on intervals 3 L), thereby eliminating its zero modes as 
well. 

As a result of the three-fold increase in the number of 
components of the quark field, unquenched simulations 
using TBC will be roughly three times as expensive as 
simulations using PBC for the same action (although this 
cost is far offset by the reduction in finite volume effects). 
To compensate for the extra fermion copies we include an 
additional factor of 1 /3 in the fermion force term (on top 
of the 1/4 to reduce the four staggered “tastes” to a single 
effective quark flavor); we therefore continue to use n/ to 
denote the total number of quark “flavors.” 

Figure |21 shows a scatter plot and run-time history of 
the temporal Polyakov for simulations done with TBC, 
using the same simulation parameters that were used 
with PBC to generate Fig. |TJ The suppression of Z( 3) 
transitions when TBC are used is striking. In fact, we 
have generated about one million configurations with 
TBC (at a larger volume V = 8 4 ), and not a single 
transition away from the “trivial” Z{ 3) phase has been 
observed. 


B. Simulation algorithms 

Dynamical simulations with staggered quarks have 
generally been done with the R algorithm. The major 
disadvantage of this algorithm is that measured quanti¬ 


ties have leading 0(At 2 ) errors ^3, where At is the step- 
size in the molecular-dynamics evolution equations. The 
step-size errors cannot be corrected by a Monte Carlo 
accept/reject procedure, because the fermion force is not 
computed explicitly in the R algorithm, rather a noisy 
estimator is used. This leads to a large change in the 
system energy during the molecular-dynamics updates, 
making the acceptance rate small if one would impose a 
Monte Carlo accept/reject step at the end of the evolu¬ 
tion. Therefore simulations are generally done at several 
values of At, and the results are extrapolated to zero 
step size. This is the method we use to simulate the 
“asqtad” actions, for which ensembles were generated at 
four values of At at each bare lattice coupling (see Ta¬ 
ble EH • The Wilson loop data at a given bare coupling 
were fit to an expansion in At, including terms up to 
0(At 6 ), excluding the linear term (constrained-curve fits 
were used, with priors for each term in the fit set to 0±5, 
see Sect. imm Some typical results for the step-size ex¬ 
trapolations are shown in Fig. [21 

Step-size errors can be eliminated using the rational 
hybrid Monte Carlo (RHMC) algorithm to handle the 
fractional powers of the staggered matrix [a! [20|. In the 
RHMC algorithm the fermion force for staggered quarks 
is computed explicitly by approximating the (jiff 4)-th 
root (or the (n//12)-th root in the case of TBC) of the 
fermion matrix with a rational function (and retaining 
only even-site couplings in M' M, as usual, to avoid a 
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FIG. 3: Representative step-size extrapolations of the fi-algorithm simulation results for the “asqtad” actions. 


further redoubling of fermion degrees-of-freedom). This 
allows an explicit computation of the force term, hence 
the step-size errors can be corrected efficiently using a 
Monte Carlo accept/reject step. 

We tested the efficiency of the RHMC algorithm by 
applying rational function approximations of various de¬ 
grees [iVratWrat]: that is, the numerator and denomina¬ 
tor are taken to be polynomials of degree N lat . This 
was done to approximate both a; -1 / 12 (for TBC with 
n/ = 1) and a; -1 / 4 (for TBC with n/ = 3), in the 
interval x £ [A m i n ,A max ], where we estimate the lower 
and upper bounds on the eigenvalue spectrum of M'M 
from the free-held values, A m ; n = (2moa) 2 and A max = 
64+ (2moa) 2 , respectively. 

Results for degrees N rat = 6, 8, 10 and 12 are shown 
in Fig. 0J for two values of the bare coupling. The max¬ 
imum error in the rational approximation, over the req¬ 
uisite eigenvalue range, falls below about 10 -5 at around 
Wat = 10, and it appears as well that at larger orders 
there is little change in the ensemble averages of the pla- 
quette, within the statistical errors. We therefore use 
Wat = 10 throughout the rest of this study (on the other 
hand, only a 5% difference in performance is observed 
for each additional order in the approximation). The 
RHMC code is found to be about two-times slower com¬ 


pared to the R algorithm, using the same At and number 
of molecular-dynamics steps. Given that the acceptance 
rate varies from about 90% at (3 = 11 to about 50% at 
(3 = 47 (see Table |1J , we conclude that the RHMC al¬ 
gorithm is about 2-4 times more expensive; however, we 
generated four sets of ensembles at different step-sizes 
for the R algorithm, in order to extrapolate to At = 0, 
hence we find that the two algorithms are, in practice, of 
comparable cost for unimproved-staggered quarks. 


C. Other algorithmic systematics 


We have explicitly computed autocorrelation times for 
the simulation data. This is important because gauge- 
field fluctuations are suppressed at weak couplings, which 
should lead to longer autocorrelation times as the cou¬ 
pling is decreased; moreover, little information on auto¬ 
correlations in the weak-coupling phase is available from 
other studies. The autocorrelation for a set of measure¬ 
ments Oi of some observable is defined, as usual, by 


Corr(Wkip) 


EiiOi - om,+N skip - o) 

Ei(Oi - d ) 2 


(18) 
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FIG. 4: Average plaquettes computed using different degrees N la ± of the polynomials in the rational approximation to the 
unimproved staggered-quark force term. Results are shown for n/ = 1 at /3 = 11 and 47. 
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FIG. 5: Autocorrelation plots for the 3x3 loops for the “asqtad” actions, at two different couplings. Configurations were 
generated by the R algorithm, with At = 0.01 and 50 molecular-dynamics steps per trajectory. 
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where N s ki p is the number of “skipped” trajectories, 
and where the autocorrelation time r is estimated from 
Corr(JV s kip) ~ e~ Nskip / T . The autocorrelation time is 
observable-dependent, and one generally expects to have 
longer autocorrelation times for a “larger” observable, 
such as a larger Wilson loop. Figure 0 shows the auto¬ 
correlation function for the 3x3 Wilson loop, the largest 
considered here, at two couplings, for the “asqtad” ac¬ 
tions. Note that the autocorrelation time is longer for 
the larger value of (3, as expected. Based on these results, 
we have chosen to skip 20 trajectories between measure¬ 
ments for the “asqtad” simulations. In the case of the 
RHMC algorithm that was used for the unimproved ac¬ 
tion, where the lowest acceptance rate is about 50%, 40 
trajectories were skipped between measurements. 

Another source of systematic error is the accuracy 
of the matrix inversion, which affects both the R al¬ 
gorithm and the RHMC algorithm. The Monte Carlo 
accept/reject step in the RHMC algorithm cannot re¬ 
move this error since matrix inversions are also required 
in computing the action for this setp. Matrix inver¬ 
sions {M'M)~ X (f) = x were done by the stabilized bi¬ 
conjugate-gradient method, with a convergence criterion 
| \{M'M)x — 4>\\ < e. We tested the precision of the in¬ 
version by comparing ensemble averages computed with 
different values of e. Results are shown in Fig. [C] for the 
average plaquette for the “asqtad” actions. Results are 
consistent within statistical errors for e < 10~ 3 ; we used 
e = ICC 5 in the production runs. 

An additional systematic error arises from propagation 
of the uncertainties in the couplings av(qh), which are 
due to the statistical errors in the measured values of 
the lxl plaquette, from which the couplings are ex¬ 
tracted. The errors in ay(gj' 1 ) propagate through to 
the couplings chv{<1*rt) that are used in the fits to the 
larger Wilson loops. This could be accounted for by in¬ 
cluding the associated uncertainty in the scale param¬ 
eters Ay (cf. (Eq. EJl), in the augmented y 2 that is 
used in the Bayesian analysis (see Sect. IIIIDI below!. 
We have instead propagated the error in the coupling 
by using a first-order approximation to the perturba¬ 
tive expansions (which is adequate for this purpose), 
— \tlWrt/^{R + T) « ci-RT otv(q* RT )- This implies that 
the statistical uncertainty A[lnWn]MC in the lxl loop 
“induces” an additional uncertainty A [In IFr^] induced 
the RxT loop, beyond its own statistical error, through 
the coupling, given by 

A [In Wrt] induced ~ fX A [In Wii]mC • (19) 

In W ii 

It turns out that the simulation error in In Wrt grows 
more rapidly than — In Wrt itself, as the loop gets larger, 
hence the “induced” error becomes less important for 
larger loops; a representative comparison of the errors 
for various loops is shown in Table El The sum of the 
two correlated errors (statistical and “induced”) is used 
in the fits for each Wilson loop. 


Loop 

-(In Wat) 

A [In Wrt] mc 
(x 10~ 5 ) 

A [111 induced 

(xl0“ 5 ) 

lxl 

0.13976 

2.2 

- 

1x2 

0.24382 

4.6 

3.8 

1x3 

0.34158 

7.9 

5.3 

2x2 

0.39270 

8.8 

6.1 

2x3 

0.52235 

13.7 

8.1 

3x3 

0.66857 

20.2 

10.4 


TABLE V: Comparison of simulation error and “induced” er¬ 
ror (due to the uncertainty in the coupling in the perturbative 
expansion), for various Wilson loops, for the unimproved ac¬ 
tions at (3 = 16 with rif = 1. 


Dependence 

on N (a = 

= 1.5) 

Cn 

JV = 4 

N = 5 

N = 6 

Cl 

1.4334(6) 

1.4334(6) 

1.4334(6) 

C2 

-1.37(4) 

-1.37(4) 

-1.37(4) 

C3 

-0.91(56) 

-0.91(56) 

-0.91(56) 

C4 

-0.1(15) 

-0.1(15) 

-0.1(15) 

C5 

- 

0.0(15) 

0.0(15) 

C6 

- 

- 

0.0(15) 

Xaug/^of 

0.56 

0.56 

0.56 



Dependence on 0 

(N = 6) 


Cn 

a = 0.5 

a = 1.0 

a = 1.5 

a = 5.0 

Cl 

1.4338(4) 

1.4335(5) 

1.4334(6) 

1.4333(6) 

C2 

-1.40(3) 

-1.38(4) 

-1.37(4) 

-1.36(4) 

C3 

-0.49(38) 

-0.79(51) 

-0.91(56) 

-0.98(77) 

C4 

0.0(5) 

-0.1(10) 

0.0(15) 

-0.5(50) 

C5 

0.0(5) 

0.0(10) 

0.0(15) 

0.0(50) 

C6 

0.0(5) 

0.0(10) 

0.0(15) 

0.0(50) 

Xaug/^of 

2.8 

0.93 

0.56 

0.28 


TABLE VI: Fit results for the 2x2 Wilson loop, for the 
“asqtad” actions with n/ = 1. The upper table shows the 
dependence of the results on the number of terms N in the 
fit, for a fixed prior width a = 1.5, while the lower table 
shows the dependence on a for fixed N = 6 . The augmented 
X 2 per degree-of-freedom (dof) is also shown. The central 
values for the prior coefficients are c„ = 0. Diagrammatic 
perturbation theory yields ci = 1.4339(0), C 2 = —1.400(2) 
and C 3 = —0.52(7) 0. More accurate Monte Carlo estimates 
of C 2 and C 3 are obtained in Sect. ITY1 using diagrammatic 
perturbation theory to constrain the lower-order terms. 


D. Fitting and truncation errors 

We use constrained curve fitting [28] in our fits to 
Eq. J5J, which allows one to incorporate the assumption 
of a convergent perturbation series in a natural way; in 
particular, a large number of higher-order terms can be 
included in the fit, without spoiling the quality of the 
results for lower-order terms that can be resolved by the 
data. 

Constrained curve fitting is motivated by Bayesian 
analysis; we refer the reader to Ref. 28] for an overview. 
In practice, we minimize an “augmented” least-squares 
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fit function Xa Ug > g iven b y 


IV. RESULTS 


+ ( 20 ) 

71=1 Gn 

where x 2 is the usual weighted sum-of-squared errors in 
fit to the Monte Carlo data, and N is the number of 
terms in the fit (cf. Eq. ©), which we take to be larger 
than the order that we anticipate can be resolved by the 
data. Fitting the data with x 2 U g favors c n ’s in the interval 
c n ± a n , which are known collectively as “priors.” 

The values of c n and a n are to be chosen based on 
theoretical expectations. In the present case one expects 
that c n = 0 ( 1 ), if the perturbative expansion is conver¬ 
gent. We therefore take c. n = 0, and use a single prior 
“width” for all orders, a n = a, which should be of 0 ( 1 ). 
The optimal value for a prior such as a can sometimes be 
determined from the data itself, by maximizing the prob¬ 
ability of obtaining the data, given the prior information, 
as a function of the prior itself; this is the so-called “em¬ 
pirical” Bayes method (for details see Ref. ,2Sj). 

The sensitivity of the data to a given order in the ex¬ 
pansion is reflected in the fit results; coefficients that are 
well determined by the data are relatively insensitive to 
changes in a and in the number of terms N in the fit, 
while fit results for coefficients that are poorly or not at 
all constrained by the data simply reproduce the priors. 
We expect that the statistical quality of the simulation 
data here should allow for useful determinations of the 
three leading orders in the perturbative expansion. 

To illustrate the quality of the constrained fits, some 
representative results are given in Table IVT1 here for the 
2x2 Wilson loop for the “asqtad” actions with n/ = 
1. The effects of varying the number of terms N in the 
fit, and of varying the prior width a, are shown. We 
see that the results for the three leading orders in the 
perturbative expansion, Ci, C 2 , and C 3 , are very insensitive 
to the details of the priors, while the fit returns the prior 
information for the fourth- and higher-order terms, which 
simply indicates that the data are not accurate enough 
to resolve those terms, as anticipated. 

The results for the three leading orders are in excel¬ 
lent agreement with diagrammatic perturbation theory, 
within the fit errors, of a few parts in 10 4 for ci, and a few 
percent for C 2 ; a reasonable NNLO signal C 3 « — 1 is also 
obtained, which is remarkable, given how little a priori 
information went into the fit (more accurate results for 
C 2 and C 3 are obtained in Sect. eh using diagrammatic 
perturbation theory to constrain the lower orders). 

The results for the Xaug °f the fits also suggest that 
an “optimal” prior width can be determined from the 
data, with a ~ 1 . 0 - 1 .5, which we also find using the 
optimization procedure given by the “empirical” Bayes 
method, alluded to above. Hence the simulation data are 
indeed consistent with the expectation that perturbation 
theory is reliable. For the rest of the fits in this paper we 
use TV = 6 , c„ = 0, and a = 1.5. 


A. Fitting strategy 


Simulation results for the perturbative coefficients of 
various Wilson loops are presented for the unimproved 
and “asqtad” actions in the next two subsections. The 
results are compared with a recent NNLO analysis using 
diagrammatic perturbation theory, which was done in the 
infinite-volume limit, and for zero quark mass, in Ref. [t|. 

We perform several types of fits for each set of actions. 
Fits are first done without input from diagrammatic per¬ 
turbation theory for a given Wilson loop (though we al¬ 
ways assume the relevant scales q RT , as well as the co¬ 
efficients for the lxl loop, in order to extract the nu¬ 
merical values of the ay couplings from the simulation 
data). We then set the first-order coefficients to their 
values from diagrammatic perturbation theory, so as to 
improve the Monte Carlo estimates of the second- and 
third-order terms, and then we similarly constrain both 
the first- and second-order coefficients, so as to obtain 
the most accurate estimate possible of the third-order 
terms. This provides us with very stringent tests of the 
Monte Carlo method, and a very precise cross-check of 
the difficult NNLO diagrammatic calculations in Ref. [?} • 
For fits that use diagrammatic values to constrain 
lower-order coefficients, we computed the leading-order 
Feynman diagram loop integrals with TBC, for the 8 4 
volume and quark masses that were used in the simula¬ 
tions, in order to consistently account for finite-volume 
and finite-mass effects in the simulation data. However 
these effects are negligible except for the largest Wil¬ 
son loops and at the smallest couplings considered here 
[the leading finite-volume corrections are ~ ay/Vol = 
0( 10 -4 ) x ay, while the leading mass-dependent correc¬ 
tions are ~ cty(moa ) 2 = 0(0.01) x a^; see also Ref. [Dj 
for an extensive finite-volume analysis in the pure-gauge 
Monte Carlo perturbation theory.] 

We also plot the following three residuals 


Ki = 


k 2 = 


— In Wrt 


a v (q* RT ) [2(R + T)\ ’ 
— In W R t 
2 (R + T) 


1,2 (Qrt) 


1 V\1RT 


CyRTOty 


( 21 ) 

( 22 ) 


and 


1 v(Qrt) 


- In W R t 2 

9(R + T) ~ Cl ’ RT<Xv ~~ c 2;RT a v 


(23) 


which provide useful visualizations of the quality of the 
data (diagrammatic perturbation theory is used for the 
coefficients in K 2 and K 3 ). We plot the residuals versus 
the appropriate coupling ay(q RT ), and in the limit of 
zero coupling one should find n n c n for the particular 
Wilson loop. Likewise a visible slope in the residual n n 
at small couplings reveals the sensitivity of the data to 
the next-order coefficient c„+ 1 . 
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FIG. 7: Plots of the simulation results for k, 1 , ^ 2 , and K 3 , for the 2x2 Wilson loop, for the unimproved actions with rif = 1. 
The simulation parameters are given in Table Q The intercepts marked by open squares (□) in the plot for each K, n give the 
results of the fits to the corresponding Cn, while the diagrammatic values [3) are indicated by filled circles (•). 


Number of flavours: rif = 1 

Monte Carlo method Diagrammatic values 

Loop 

Cl 

C 2 

C3 

Cl 

C 2 

C3 

1 x 2 

1x3 

2 x 2 

2x3 

3x3 

1.2038(4) 

1.2586(5) 

1.4334(6) 

1.5163(7) 

1.6080(8) 

-1.327(29) 

-1.253(36) 

-1.368(39) 

-1.281(47) 

-1.230(56) 

-1.26(46) 

-1.43(56) 

-0.95(59) 

-1.11(71) 

-0.45(82) 

1.2039(0) 

1.2589(0) 

1.4339(0) 

1.5172(0) 

1.6090(0) 

-1.335(0) 

-1.277(1) 

-1.400(2) 

-1.351(3) 

-1.298(6) 

-1.10(3) 

-0.95(6) 

-0.52(8) 

-0.26(12) 

0.67(25) 


Number of flavours: rif = 3 

Monte Carlo method Diagrammatic values 

Loop 

Cl 

C 2 

C3 

Cl 

C 2 

C3 

1 x 2 

1x3 

2 x 2 

2x3 

3x3 

1.2039(4) 

1.2590(5) 

1.4337(5) 

1.5169(6) 

1.6088(8) 

-1.480(28) 

-1.444(34) 

-1.529(38) 

-1.476(46) 

-1.431(53) 

-0.28(45) 

-0.06(54) 

0.02(60) 

0.20(71) 

0.65(81) 

1.2039(0) 

1.2589(0) 

1.4339(0) 

1.5172(0) 

1.6090(0) 

-1.485(0) 

-1.437(1) 

-1.551(2) 

-1.513(4) 

-1.463(11) 

-0.11(9) 

0 . 02 ( 11 ) 

0.51(13) 

0.73(16) 

1.62(30) 


TABLE VII: Perturbative coefficients of various small Wilson loops for the unimproved actions, with rif = 1, and with rif = 3. 
The coefficients from diagrammatic perturbation theory Q are compared with the results obtained from the Monte Carlo 
simulations. These fits to the Monte Carlo data used no diagrammatic input, except in the determination of the couplings for 
the expansion Eq. 


Diagrammatic input: ci 
rif — 1 rif = 3 

Loop 

C 2 

C3 

C 2 

C3 

1 x 2 

1x3 

2 x 2 

2x3 

3x3 

-1.332(9) 

-1.270(11) 

-1.403(12) 

-1.342(14) 

-1.292(17) 

-1.19(22) 

-1.20(27) 

-0.48(29) 

-0.28(34) 

0.38(41) 

-1.480(10) 

-1.434(12) 

-1.541(13) 

-1.499(16) 

-1.442(20) 

-0.28(25) 

-0.19(30) 

0.20(33) 

0.52(39) 

0.80(46) 


Diagrammatic input: ci ,2 
rif = 1 rif = 3 

Loop 

C3 

C3 

1 x 2 

1x3 

2 x 2 

2x3 

3x3 

-0.11(9) 

0 . 02 ( 11 ) 

0.51(13) 

0.73(16) 

1.62(30) 

-0.17(10) 

- 0 . 12 ( 11 ) 

0.42(13) 

0.81(16) 

1.16(27) 


TABLE VIII: Monte Carlo results for perturbative coefficients for the unimproved actions, where in the left table the first-order 
term is set to its diagrammatic value, while in the right table both ci and C 2 are set to their diagrammatic values. 


B. Unimproved actions 

The unimproved actions were simulated for n/ = 1 
and rif = 3. The sensitivity of the data to three lead¬ 
ing orders in the perturbative expansions is apparent in 
the plots of the three residuals «i, « 2 , and K 3 , which are 
shown for a representative Wilson loop in Fig. 0 The 


slope and curvature in the plot for K \, for instance, in¬ 
dicate the sensitivity to C 2 and C 3 , while the fact that 
the plot for K 3 has no noticeable slope indicates that the 
data are insensitive to C 4 , within the statistical errors. 

Fit results for ci, C 2 , and C 3 , without input from dia¬ 
grammatic perturbation theory (except for the couplings) 
are given in Table IvTTl for Wilson loops up to 3x3. As 
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FIG. 8: Plots of k, 1, ^2, and K 3 for the 2x2 loop with the “asqtad” actions. Simulation parameters are given in Table ITTTl 
The results of fits to the Monte Carlo data for the perturbative coefficients are indicated by open squares (□), while the 
diagrammatic results are shown by filled circles (•). 


Diagrammatic input: ci,2 

Loop 

C3 

1x2 

0.21(4) 

1x3 

0.28(5) 

2x2 

0.61(6) 

2x3 

0.84(8) 

3x3 

0.84(10) 



Monte Carlo method 

Diagrammatic values 


Diagrammatic input: ci 

Loop 

Cl 

C2 

C3 

Cl 

C2 

C3 


Loop 

C2 

C3 

1x2 

0.9251(3) 

-0.644(13) 

0.20(18) 

0.9252(0) 

-0.646(0) 

0.23(5) 


1x2 

-0.649(5) 

0.26(12) 

1x3 

0.9845(3) 

-0.599(14) 

0.37(19) 

0.9845(0) 

-0.595(1) 

0.38(6) 


1x3 

-0.600(6) 

0.39(13) 

2x2 

1.1499(4) 

-0.641(15) 

0.58(20) 

1.1499(0) 

-0.643(2) 

0.59(9) 


2x2 

-0.642(7) 

0.59(14) 

2x3 

1.2342(4) 

-0.599(19) 

0.88(26) 

1.2341(0) 

-0.595(3) 

0.85(16) 


2x3 

-0.594(7) 

0.82(15) 

3x3 

1.3235(5) 

-0.545(19) 

1.16(23) 

1.3235(0) 

-0.522(4) 

0.96(19) 


3x3 

-0.546(8) 

1.17(17) 


TABLE IX: Perturbative coefficients for various small Wilson loops for the “asqtad” actions. There is no diagrammatic input 
in the Monte Carlo results in the left-most table (except for the couplings). In the middle table the Monte Carlo results are 
shown with the first-order term set to its diagrammatic value, while in the right-most table both ci and C2 are set to their 
diagrammatic values. 


discussed in Sect. ITTT PI excellent agreement with dia¬ 
grammatic calculations is obtained, within the fit errors, 
which for Ci are typically a few parts in 10 4 , and for C 2 are 
typically a few percent. The third order terms are also re¬ 
solved, which is remarkable given that no diagrammatic 
input for these Wilson loops was used. The simulations 
also clearly resolve the perturbative dependence on the 
number of flavors, which in the case of C 2 shows a change 
of 3-4 standard deviations from ny = lton/=3. 

The accuracy of the Monte Carlo results for the higher- 
order coefficients can be dramatically improved by fur¬ 
ther constraining the lower-order terms using available 
diagrammatic input. Results are shown in Table Enn 
where we first fix C\ to its diagrammatic value, and then 
fix both Ci and C 2 ; in the first case, the fit errors in C 2 and 
C 3 are reduced by factors of about 2-3 compared to the 
results in Table lvT!l while in the second case the errors 
in C 3 are reduced by a factor of about 5. Agreement with 
the diagrammatic values is obtained in all cases, within 
these greatly reduced errors, and in most cases the Monte 
Carlo results have errors that are comparable to those in 
the diagrammatic evaluations. 


C. The “asqtad” actions 

Our results for the perturbative expansions of Wilson 
loops for the “asqtad” actions are shown in Fig. |B1 and 
in Table Hxl The data is sensitive to the three leading 


orders of the perturbative expansion, and the agreement 
with diagrammatic results is again impressive. 

We note that the uncertainties in the coefficients ob¬ 
tained from these simulations are about a factor of 3 
smaller than for the unimproved results presented in 
the preceding subsection, which owes to the fact that 
the “asqtad” simulations were done at larger couplings 
[the largest coupling for the unimproved simulations was 
av(qh) ~ 0.06, while for the “asqtad” simulations it was 
a vil 11 ) ~ 0.13, cf. Tables HI and IIIII . This fact should 
clearly be heeded in future studies, where one might work 
at still larger couplings, as long as the theory remains in 
the perturbative phase, as judged for example by simula¬ 
tion measurements of the Polyakov line, which provides 
an order parameter for confinement (a systematic ap¬ 
proach to optimizing the choice of couplings is given in 
Refs. IJ, 223). 


V. CONCLUSIONS AND OUTLOOK 

Perturbative coefficients of Wilson loops were ex¬ 
tracted from unquenched QCD simulations at weak cou¬ 
plings. Two sets of actions were analyzed: unimproved 
gluon and staggered-quark actions in one set, and 0 (a 2 ) 
improved actions in the other. Simulations were also 
done for different numbers of dynamical fermions. An 
extensive analysis of systematic uncertainties was made; 
constrained-curve fitting in particular was used to extract 
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as much information as possible from the simulation data. 

The Monte Carlo results for the perturbative coeffi¬ 
cients were found to be in excellent agreement with calcu¬ 
lations using diagrammatic perturbation theory, through 
next-to-next-to leading order. Results were obtained for 
the first-order coefficients with uncertainties of a few 
parts in 10 4 , while the second-order coefficients were ob¬ 
tained to a few-percent precision, without any input from 
diagrammatic perturbation theory (except for the pertur¬ 
bative expansion of the lxl plaquette, which was used 
to extract the relevant couplings ay(q*) from the simu¬ 
lation data). The results also show that the Monte Carlo 
perturbation theory is clearly sensitive to the number 
of dynamical fermion. Furthermore, when the two lead¬ 
ing perturbative coefficients were constrained to their di¬ 
agrammatic values, the third-order term was obtained 
with a precision comparable to that from the NNLO dia¬ 
grammatic analysis |2j, which required the evaluation of 
about fifty Feynman diagrams. 

These results provide a stringent test of the Monte 
Carlo method as applied to highly-improved lattice ac¬ 
tions with dynamical fermions, and also provide an im¬ 
portant high-precision cross-check of the perturbation 
theory input to a recent determination of ajj^(Mz) by 
the HPQCD collaboration j§|. 

Remarkably little computational power was needed, 
since the simulations were done on small volumes ( 8 4 ), 
thanks to the use of twisted boundary conditions, which 


eliminate lattice zero modes, and which suppress other 
finite-volume effects; the entire set of simulations for the 
“asqtad” actions only required the equivalent of about 60 
months of run-time on a single 3 GHz processor. Larger 
lattices might be needed for other quantities, such as 
quark propagators Q that are not as dominated by ul¬ 
traviolet modes as the small Wilson loops anal yzed here, 
though an earlier study for pure-gauge theories |l3| found 
that finite-volume effects in such quantities could be re¬ 
moved by working on several relatively small volumes. 

Additional work using this method is in progress, in¬ 
cluding the determination of the NNLO mass renormal¬ 
ization for the NRQCD heavy-quark action, and prelimi¬ 
nary investigations of currents with infrared singularities. 


Acknowledgments 

We thank Peter Lepage, Christine Davies, and 
Matthew Nobes for fruitful conversations. We also thank 
Peter Lepage for providing us with his Bayesian-analysis 
software. We are grateful to Randy Lewis for providing 
resources including the use of his computer cluster at the 
University of Regina. Computations were also performed 
on facilities provided by WestGrid l32i| . This work was 
supported in part by the Natural Sciences and Engineer¬ 
ing Research Council of Canada. 


[1] C. T. H. Davies et al, Phys. Rev. D 56, 2755 (1997). 

[2] Q. Mason et al, HPQCD Collaboration, Phys. Rev. Lett. 
95, 052002 (2005). 

[3] H. D. Trottier, Nucl. Phys. Proc. Suppl. 129, 142 (2004). 

[4] Q. J. Mason, Ph.D. thesis, Cornell University (2004). 

[5] M. Liischer and P. Weisz, Nucl. Phys. B 266,309 (1986). 

[6] M. Liischer and P. Weisz, Nucl. Phys. B 452, 234 (1995); 
B. Alles, M. Campostrini, A. Feo and H. Panagopou- 
los, Phys. Lett. B 324, 433 (1994); B. Alles, A. Feo and 
H. Panagopoulos, Nucl. Phys. B 491, 498 (1997); E. Fol- 
lana and H. Panagopoulos, Phys. Rev. D 63, 017501 
( 2001 ). 

[7] Q. Mason and H. D. Trottier, in preparation. 

[8] Q. Mason, H. D. Trottier, R. Horgan, C. T. H. Davies 
and G. P. Lepage, hep-ph/0511160 

[9] G. P. Lepage, Phys. Rev. D 59, 074502 (1999). 

[10] K. Symanzik, Nucl. Phys. B 226, 187 (1983). 

[11] K. Orginos and D. Toussaint, Phys. Rev. D 59, 014501 
(1999); C. W. Bernard et al, MILC Collaboration, Phys. 
Rev. D 64, 054506 (2001). 

[12] W. C. Dimm, G. P. Lepage and P. B. Mackenzie, Nucl. 
Phys. Proc. Suppl. 42, 403 (1995). 

[13] H. D. Trottier, N. H. Shakespeare, G. P. Lepage and 
P. B. Mackenzie, Phys. Rev. D 65, 094502 (2002). 

[14] K. J. Juge, Nucl. Phys. Proc. Suppl. 94, 584 (2001); 106, 
847 (2002). 

[15] A. Hart, R. R. Horgan and L. C. Storoni, Phys. Rev. D 
70, 034501 (2004). 

[16] Part of this work was presented in K. Y. Wong, Ph.D. 


thesis, Simon Fraser University, 2005. 

[17] F. Di Renzo et al, Nucl. Phys. B 426, 675 (1994). 

[18] F. Di Renzo and L. Scorzato, JHEP 0410, 073, (2004). 

[19] G. ’tHooft, Nucl. Phys. B 153, 141 (1979). 

[20] G. Parisi, LNF-84/4-P, Invited talk given at Summer 
Inst. Progress in Gauge Field Theory, Cargese, France, 
Sep 1-15 (1983). 

[21] R. Wohlert, DESY 87/069 (1987). 

[22] A. Gonzalez Arroyo and C. P. Korthals Altes, Nucl. Phys. 
B 311, 433 (1988). 

[23] G. P. Lepage and P. B. Mackenzie, Phys. Rev. D 48, 2250 
(1993). 

[24] S. J. Brodsky, G. P. Lepage and P. B. Mackenzie, Phys. 
Rev. D 28, 228 (1983). 

[25] I. Horvath, A. D. Kennedy and S. Sint, Nucl. Phys. Proc. 
Suppl. 73, 834 (1999). 

[26] M. A. Clark and A. D. Kennedy, Nucl. Phys. Proc. Suppl. 
129, 850 (2004). 

[27] S. Gottlieb et al, Phys. Rev. D 35, 2531 (1987). 

[28] G. P. Lepage et al, Nucl. Phys. Proc. Suppl. 106, 12 

( 2002 ). 

[29] Y. Schroder, Phys. Lett. B 447, 321 (1999). 

[30] We recomputed the leading-order Feynman diagram loop 
integrals here for finite-volume (8 4 ) TBC perturbation 
theory, and for the finite quark mass; see Sect. EH 

[31] M. G. Alford, W. Dimm, G. P. Lepage, G. Hockney and 
P. B. Mackenzie, Phys. Lett. B 361, 87 (1995). 

[32] See http: //www. westgrid. ca/home . html 





